Phylogeography and molecular diversity of two highly abundant Themisto amphipod species in a rapidly changing Arctic Ocean

Abstract Rapid warming in the Arctic is drastically impacting marine ecosystems, affecting species communities and food‐web structure. Pelagic Themisto amphipods are a major component of the Arctic zooplankton community and represent a key link between secondary producers and marine vertebrates at higher trophic levels. Two co‐existing species dominate in the region: the larger Themisto libellula, considered a true polar species and the smaller Themisto abyssorum, a sub‐Arctic, boreal‐Atlantic species. Recent changes in abundance and distribution ranges have been detected in both species, likely due to the Atlantification of the Arctic. The ecology and genetic structure of these species are understudied, despite their high biomass and importance in the food web. For both species, we assessed genetic diversity, patterns of spatial genetic structure and demographic history using samples from the Greenland shelf, Fram Strait and Svalbard. This was achieved by analysing variation on the mitochondrial cytochrome c oxidase subunit 1 gene (mtCOI). The results revealed contrasting levels of mtCOI diversity: low levels in T. libellula and high levels in T. abyssorum. A lack of spatial genetic structure and a high degree of genetic connectivity were detected in both species in the study region. These patterns of diversity are potentially linked to the impacts of the Last Glacial Maximum. T. libellula populations may have been isolated in glacial refugia, undergoing gene flow restriction and vicariant effects, followed by a population expansion after deglaciation. Whereas T. abyssorum likely maintained a stable, widely distributed metapopulation further south, explaining the high diversity and connectivity. This study provides new data on the phylogeography of two ecologically important species, which can contribute to predicting how zooplankton communities and food‐web structure will manifest in the rapidly changing Arctic.


| INTRODUC TI ON
The Arctic is currently warming three to four times faster than the global mean, and many consequences of climate change have now manifested throughout the region (Pörtner et al., 2019;Rantanen et al., 2022). Ocean and air temperatures are increasing rapidly, sea ice coverage and thickness are declining, and glaciers are melting region wide (Stroeve et al., 2012;Wang et al., 2020). These environmental changes are already having drastic impacts on the marine ecosystem, affecting species composition, distribution and food-web structure in the Arctic Ocean (Gluchowska et al., 2017;Weydmann et al., 2014).
The inflow of increasingly warmer Atlantic water into the high Arctic is a major driver of the phenomenon commonly referred to as the 'Atlantification' of the Arctic (Polyakov et al., 2017). The Fram Strait, between Greenland and the Svalbard Archipelago, is known as 'the gateway to the Arctic' and is the largest source of oceanic heat into the Arctic Basin (Beszczynska-Moeller et al., 2011). The heat exchange between water masses in the Fram Strait is influenced by two opposing currents: the north-bound West Spitsbergen Current (WSC) and the south-bound East Greenland Current (EGC) (Figure 1). The WSC carries warmer and more saline Atlantic water northwards into the Arctic Basin, while the ECG carries cold, fresher water and sea ice southwards through the Fram Strait ( Figure 1) (Wang et al., 2019).
The WSC transports increasing amounts of nutrients and subarctic and boreal-Atlantic plankton species into the Arctic (Csapó et al., 2021;Descôteaux et al., 2021;Gluchowska et al., 2017;Kraft et al., 2013), the consequences of which are not yet widely studied. Zooplankton community composition and functional changes have been detected in the region and linked to this Atlantification (Gluchowska et al., 2017). Poleward range expansions of subarctic and boreal-Atlantic zooplankton species, as well as poleward contractions of polar species have also been observed (Andrews et al., 2019;Basedow et al., 2018;Csapó et al., 2021;Dalpadado et al., 2016).
Environmental changes will continue to alter the distribution of suitable habitats and are expected to have an impact on gene flow and genetic structure in a wide range of zooplankton species in the Arctic and its marginal seas (Hardy et al., 2011;Tempestini et al., 2020).
Pelagic amphipods, and Hyperiidea in particular, are dominant planktonic crustaceans in terms of biomass in the polar regions and are a crucial component of the Arctic food web (Vinogradov et al., 1996). They are an important link between secondary producers (mesoplanktonic grazers) and marine vertebrates at higher trophic levels. Their role as prey in the Arctic zooplankton community has been described as on par with that of krill and copepods (Bowman, 1960;Dalpadado, 2002). considered a genuine polar species and Themisto abyssorum (Boeck, 1871), a subarctic, boreal-Atlantic species. An invasive third species from the North Atlantic, Themisto compressa (Goes, 1865), has also been reported in the Fram Strait, albeit in lower abundances (Kraft et al., 2013). Both T. libellula and T. abyssorum are preyed upon by fish, seabirds and marine mammals throughout the Arctic (reviewed in . They are visual predators of mesoand microzooplankton and although their geographical distributions overlap, they occupy different ecological niches (Auel et al., 2002).
Themisto libellula (Figure 2a) is a cold-adapted species and its distributional range includes the Central Arctic Basin as well as its marginal seas . It can grow up to 60 mm in size and has a life cycle of up to 4 years (Auel & Werner, 2003;Kraft, 2010). An important component of its diet are ice-dependent, herbivorous copepods, indicating a reliance on the cryo-pelagic pathway (Auel & Werner, 2003;Kohlbach et al., 2016). It is a key prey item for seabirds, including the little auk (Alle alle (Linnaeus, 1758)), Arctic fish species such as polar cod (Boreogadus saida (Lepechin, 1774)) and commercially important fish such as Atlantic cod (Gadus morhua Linnaeus, 1758) and salmonid species (reviewed in Pinchuk et al., 2013).
Themisto abyssorum (Figure 2b) is largely found in waters of Atlantic origin (Mumm et al., 1998) in the marginal Arctic seas and has a life cycle of up to 2 years . It is smaller than T. libellula, and while lipid content by body mass is comparable, its smaller body size makes it an overall less nutritious prey (Auel et al., 2002). Previous studies indicate that the diet of T.
abyssorum is less sea ice-dependent than that of T. libellula, with a wider prey spectrum and a higher trophic position in the Arctic food web (Auel et al., 2002;Kohlbach et al., 2016).
Both species have exhibited recent changes in abundance and distribution, likely as a result of the Atlantification of the Arctic.
Themisto libellula has decreased in abundances, whereas T. abyssorum has become more abundant in the Fram Strait and the Barents Sea (CAFF, 2017;. These changes in abundance and distribution could have strong implications at higher trophic levels, with the loss of the highly nutritious T. libellula as a key prey item for many species (reviewed in . Many aspects of the current genetic diversity, population structure and phylogeography of T. libellula and T. abyssorum are not well studied, despite their importance in the Arctic food web and biogeochemical cycles . Investigating patterns of genetic diversity over a spatial distribution can give insight into how historic evolutionary events have impacted genetic variation in a species, as well as ongoing processes such as gene flow (Nei, 1987;Slatkin, 1987). The phylogenetic relationships between different Themisto species have been addressed (Tempestini et al., 2017), but our understanding of genetic diversity and distribution patterns are still limited within species, due to the lack of studies on small-scale geographical structure. This study is the first, to our knowledge, to focus in depth on the genetic diversity and structure of T. abyssorum, and the first to compare these two congeneric species using molecular barcoding.
We applied the commonly used 'barcode' gene for animal taxa, a fragment of the mitochondrial cytochrome c oxidase subunit 1 gene (hereafter mtCOI) (Hebert et al., 2003), to characterize the phylogeography and molecular diversity of the two congeners T. abyssorum and T. libellula. We aimed to (i) assess and compare the levels of intra-and interspecific genetic diversity, (ii) determine present spatial patterns of genetic structure and connectivity, as well as (iii) provide an overview of the demographic history of the two Themisto species, across a broad geographical gradient in the Atlantic sector of the Arctic. This was achieved by analysing and comparing mtCOI gene variation according to the geographical location of samples from the Greenland shelf, the Fram Strait and the Svalbard Archipelago.
Considering their relationship as congeneric species (Tempestini et al., 2017), similar life history patterns and dispersal capacities, we hypothesized comparative levels of genetic diversity at the mtCOI region. Both species have a holoplanktonic lifestyle, and considering F I G U R E 2 (a) Themisto libellula and (b) Themisto abyssorum. Size of the adult individuals typically ranged from 0.6 to 2.5 cm for T. libellula and from 0.4 to 1.3 cm for T. abyssorum in the study area (C. Havermans, personal observation).
the multidirectional currents and the lack of major geographical barriers in the study region, we predicted relatively high levels of genetic connectivity. As an alternative hypothesis, divergent genetic diversity patterns could be expected between the two species based on their current distributions and their link to distinct water masses (Dalpadado, 2002). From an eco-evolutionary perspective, both species may have undergone different gene flow restrictions during the last glacial maximum (LGM) (Hardy et al., 2011). Themisto libellula may have had to survive the ice-covered glacial periods in Arctic refugia, leading to diversification events after a bottleneck, whereas T. abyssorum may have been pushed south out of the Arctic, surviving as one large, connected population in ice-free waters. Hence, alternatively, we expected potentially high levels of genetic divergence, with the possibility of cryptic species within T. libellula, and a higher genetic homogeneity in T. abyssorum.

| Sample collection
Zooplankton samples containing Themisto amphipods were collected during oceanographic research cruises between 2016 and 2020 (Table 1 and Figure 1). Specimens from the Fram Strait and East Greenland were collected on the R/V Polarstern (Knust, 2017) cruises PS100 in 2016 (Kanzow, 2017) and PS107 in 2017 (Schewe, 2018), with a 150 μm Multinet and 300 and 500 μm Bongo nets. Multinet hauls were carried out vertically in the water column (0.5 m/s), Bongo nets were towed obliquely at 2 knots ship's speed, with a wire length varying between 20 and 450 m. Additional Fram Strait and East Greenland specimens were collected during the cruise TUNU-VII in September 2017 on the R/V Helmer Hanssen, with a Campelen 1800 shrimp trawl (Walsh & McCallum, 1997). Specimens from the Svalbard Archipelago were collected in 2020 on the R/V Heincke cruise HE560 (Mark et al., 2022), using 300 μm and 500 μm Bongo nets and a pelagic trawl net fitted with a fish lift (Holst & McDonald, 2000). Samples were preserved immediately after collecting and sorting procedures in 96% undenatured ethanol. Six major geographical regions were defined based on ocean bathymetry and influences of major oceanic currents (Figure 2). The target gene fragment of this study was the mtCOI barcoding region. This gene fragment has the advantage that it has low or no recombination, maternal inheritance and a faster evolutionary rate compared to nuclear DNA (Moritz et al., 1987). Diversity on mtCOI is higher within species than between species, making it a suitable tool for estimating intraspecific diversity in congeneric species (Meyer & Paulay, 2005;Vieira et al., 2022). Furthermore, population bottlenecks and expansions leave signatures in the mtCOI region, making it a popular tool for inferring elements of a species demographic history (Beermann et al., 2020;Hewitt, 2000). The 658 base pair (bp) barcoding region of the mtCOI was amplified using the universal 'Folmer' primers HCO 2198 (5′-TAAAC TTC AGG GTG ACC AAA AAATCA-3′) and LCO 1490 (5′-GGTCAACAAAT-CTAAAGATATTGG -3′) (Folmer et al., 1994). The reaction mixture had a total volume of 25 μL and consisted of: 1x PCR buffer, 0.2 mM dNTPs, 0.5 μM forward and reverse primers, 0.2 U/μL HotMaster Taq DNA polymerase (QuantaBio), 1 μL (ca. 30 ng) of DNA template and molecular-grade water. DNA template was substituted with molecular-grade water in negative controls. The PCR amplification was carried out using the following programme: initial denaturation at 95°C for 2 min; followed by 36 cycles of denaturation at 94°C for 20s, annealing at 42°C for 20s, extension at 65°C plus a final extension at 65°C for 15 min. PCR products were checked for quality and length using electrophoresis on GelRed-stained, 2% agarose gel and were bidirectionally sequenced by EUROFINS (Germany).  Note: Numbers in parentheses () indicate shorter sequences (<658 bp) that were included in the genetic distance analysis and ML tree but removed from subsequent analysis.

| Genetic diversity analysis and demographic history
Sequences were collapsed into unique haplotypes and used to construct a maximum likelihood (ML) phylogeny tree using PhyML was used to edit the subsequent ML tree. Statistical support for the clades was estimated with 1000 bootstrap replicates.
The T. compressa sequences (N = 22) were only included in the genetic distance analysis and the ML tree to illustrate the relationship between the three congeneric species but were not used further in this study. An additional 33 publicly available sequences of T.
libellula from the Pacific sector of the Arctic (Tempestini et al., 2020) (GenBank accession numbers: MT832367-MT832393) were included in the genetic distance analysis and ML tree to investigate their relationship to the T. libellula in the present study, but not included in subsequent analysis. Furthermore, five sequences of the hyperiid Hyperiella dilatata Stebbing, 1888 from a previous study  were added to the tree as an out-  (Rozas et al., 2017). The relationships between identified haplotypes were explored through haplotype networks created in PopART 1.7 software (Leigh & Bryant, 2015), using the Templeton, Crandall and Sing (TCS) method (Clement et al., 2002).
The TCS method is based on a maximum parsimony algorithm.
Three different classes of common analyses were used to infer demographic history. Tajima's D (Class I), which is based on the frequency of segregating nucleotide sites (Tajima, 1989) and Ramos-Onsins and Rozas's (R 2 , Class I), which is calculated using the number difference between singleton count and the average number of nucleotide differences (Ramos-Onsins & Rozas, 2002).
Fu's F (Class II) (Fu, 1997), which uses the distribution of haplotype frequencies and, along with R 2 , is considered to be more sensitive than Tajima's D (Ramos-Onsins & Rozas, 2002). In addition, mismatch distribution analysis based on pairwise sequence differences, and Harpending's raggedness index (r, Class III) (Ramos-Onsins & Rozas, 2002) were calculated using a population growth-decline model. All demographic tests were calculated in DnaSP 6 (Rozas et al., 2017) and their significance assessed using 10,000 bootstrap replicates.

| Population structure and connectivity analysis
To test for hierarchical population genetic differentiation, an analysis of molecular variance (AMOVA) was performed using the distribution of variation at the regional, sampling site and individual levels (Excoffier et al., 1992). Genetic connectivity between geographical regions was estimated by calculating pairwise θ ST (fixation index among populations), for which a significance level of 0.05 was determined using 10,100 permutations (Holsinger & Weir, 2009). These analyses were conducted using Arlequin 3.5.2.

| Maximum likelihood clustering of T. abyssorum and T. libellula
All Themisto barcodes clustered in the ML tree with high bootstrap support for species-level clusters (>85%) (Figure 3)
Intraspecific genetic distance was the lowest among the species studied (0.07%; Table 2). In contrast, 136 polymorphic sites and 68 parsimony informative sites were identified for T. abyssorum, as well as 115 unique haplotypes belonging to 152 individuals. Genetic diversity indices were higher with H d = 0.975 (SD = 0.008), π = 0.01218 (SD = 0.00089) and K = 8.02 (Table 3). Intraspecific genetic distance was higher than the other species clusters at 1.29% (Table 2).
TCS haplotype networks based on mtCOI region illustrate contrasting diversity between the two species in the Atlantic Arctic study region (Figure 4). T. libellula showed two dominant haplotypes and 13 singletons. All haplotypes diverge by only one or two of mutational steps, which is reflected in the low nucleotide diversity value (π; Table 3). In contrast, T. abyssorum exhibited a large number of highly divergent haplotypes, with only one haplotype occurring in more than 10 individuals, and 107 singleton haplotypes.
The results of the Tajima's D and Fu's F tests for neutrality were significantly negative for both species (Table 3)

| Population structure and connectivity
The data revealed a lack of genetic structure for both species in all the sampled geographical regions. The TCS haplotype networks for T.
libellula showed that the two most common haplotypes (Hap 1 = 197 and Hap 2 = 35 individuals; Figure 5) were present in all six of the main geographical regions sampled (Figure 4). The third most common haplotype (Hap 9 = 5) was present at three out of the six regions ( Figure 5).
This lack of spatial genetic structure was further underlined by a single haplotype (Hap 1) dominating the haplotype frequency at each sampling station, except one, for T. libellula ( Figure 5). The haplotype network for T.abyssorum showed that the dominant haplotype (Hap 8 = 23) was present at four out of five of the regions sampled (no T. abyssorum sequences were obtained from South Spitsbergen) and is present at all except the two Nordaustlandet stations. Every station, except for one Fram Strait station is dominated by singletons ( Figure 5).

TA B L E 2
Estimates of inter-and intraspecific genetic distances (above and in the diagonal respectively).
any of the variance components (among regions, among stations or among individuals). Based on these results, neither species exhibited significant geographical structure among the sampled populations (Table S1). Pairwise θ ST comparisons between the geographical regions showed no significant differences between pairs of regions, further indicating a lack of genetic structure and high levels of connectivity between the Atlantic and Arctic regions in both T. libellula and T. abyssorum (Table S2).

| DISCUSS ION
Two abundant pelagic amphipod species were analysed for genetic diversity, demographic history, and patterns of spatial genetic structure and connectivity, using variation at the mitochondrial COI region. Overall, the levels of mtCOI diversity were strikingly differ-

| Genetic diversity within T. libellula and T. abyssorum
We should be regarded as a species complex. The latter is supported by our K2P divergence estimates and has previously been suggested by Tempestini et al. (2017). Whether these species-level lineages also differ in morphology cannot be confirmed with the present data and should be the focus of future studies.
Haplotypic diversity in T. abyssorum was more than two-fold higher, and the number of unique haplotypes was close to 10 times that found in T. libellula, despite a smaller sample size. Nucleotide diversity was also more than one magnitude higher than that of T. libellula, alluding to high genetic distances between haplotypes. A higher haplotype diversity was also observed in temperate Pseudocalanus copepods when compared to an Arctic congeneric species (Questel et al., 2016). The high genetic diversity in T. abyssorum is also evident in the haplotype network, which is a 'diffuse' shape with numerous mutations separating the many haplotypes. This network shape has been found in other species with deep-sea distributions, including the benthic shrimp, Nematocarcinus lanceopes (Raupach et al., 2010) and the squat lobsters, Munida endeaviurae and Munida gracilis (Yan et al., 2020). or species may have migrated southwards to more ice-free waters. This extreme reduction of available habitat, combined with periods of isolation in refugia, caused declines in population size and genetic diversity across many marine Arctic taxa (Hardy et al., 2011;Hewitt, 2004).
Subsequent deglaciation was followed by rapid recolonization, causing founder effects in species of which few individuals, or individuals with low genetic variation, were able to recolonize the area (Hardy et al., 2011;Shimizu et al., 2018). The signatures of the last glacial cycle of the contemporary genetic diversity thus depend on the size and number of refugia, the effective population size and the standing genetic variation of the survivors (Hardy et al., 2011;Hewitt, 2000). Generally, similar genetic signatures of the LGM have been identified in other zooplankton taxa in the Arctic including the pteropod L. helicina (Shimizu et al., 2018;Sromek et al., 2015), the copepod Pseudocalanus moultoni (Aarbakke et al., 2014) as well as in fish species, such as the circum-Arctic capelin (Mallotus villosus; Dodson et al., 2007). However it remains unclear whether the low genetic diversity in certain Arctic species is caused entirely by limited former refugia, colonization following deglaciation and consequent founder effects (Shimizu et al., 2018). Mutations that have occurred since may explain the presence of rarer haplotypes (Lasota et al., 2016), however, the time since the LGM is likely insufficient for many mutations to have accumulated in T. libellula (Grant et al., 2011). An alternative explanation that cannot be eliminated based on the current dataset is that of a mitochondrial selective sweep reducing genetic variability on the mtCOI gene (e.g. Hurtado et al., 2004;Jenkins et al., 2018;Lasota et al., 2016;Strasser & Barber, 2009). This hypothesis could however be ruled out when carrying out a multigene study, as it would be unlikely that selective sweeps occur in several unlinked mitochondrial and nuclear markers (Strasser & Barber, 2009).

| Population structure and connectivity
The results herein did not reveal any significant genetic differentiation among sampling locations for either of the two species, as clearly illustrated in the haplotype networks, haplotype distribution maps and subsequent statistical analysis (AMOVA, θ ST ).
Dominant haplotypes were present across the sampling area for both species and all pairwise θ ST values were low and insignificant between the major geographical regions. Collectively, we interpret this as an indication of high levels of connectivity and genetic homogeneity for both species, within the Atlantic-Arctic region.
Weak or non-existent population structure is common among marine species with pelagic life stages, such as Themisto amphipods, due to the high potential for dispersal and gene flow through ocean currents (Hardy et al., 2011). Similar high levels of population connectivity have been demonstrated in other key Arctic zooplankton species including Calanus copepods (Weydmann et al., 2014, 2018), chaetognaths (DeHart et al., 2020, as well as in intertidal amphipods (Grabowski et al., 2019), and for multiple polychaete and echinoderm species with planktonic life stages (Hardy et al., 2011). Increases in presence and abundances of T. libellula were observed in colder years in the southern Bering Sea, and found to disappear in warmer years (Pinchuk et al., 2013). Mass occurrences of T. libellula have been reported in northern parts of the Bering Sea during periods of cold-water inflow, but decreased with subsequent warming (Volkov, 2012). Although these distribution changes are not all poleward, the authors identify temperature changes as the main driver behind species presence and abundance, with T. libellula populations closely following the movements of colder waters.  (Kraft et al., 2013;Schröter et al., 2019). Other Arctic pelagic zooplankton species, including krill species (Buchholz et al., 2010) and Calanus copepods (Weydmann et al., 2014), have shown community and distribution changes as a result of Atlantification around the Svalbard Archipelago and in Fram Strait. These changes in distribution and abundance associated with the changing environmental conditions in the Arctic are evidence that zooplankton communities are already in transition to a warmer Arctic dominated by boreal species (Csapó et al., 2021).
A combination of factors, including low genetic diversity at the COI gene, long life cycle , cold-water affinity and its potential associated adaptions (Pinchuk et al., 2013), and reliance on the cryo-pelagic pathway (Auel et al., 2002;Kohlbach et al., 2016), may well lead to T. libellula being negatively impacted by the Atlantification of the Arctic. This could result in loss of intraspecific diversity, local extinctions and further poleward distribution contractions with unknown effects on the food web. In contrast, T.
abyssorum could benefit from an increasing Atlantification, as it displays an Atlantic affinity and a shorter life cycle . Moreover, the high standing genetic variation observed in T. abyssorum, if confirmed at the genome level, can be an important source of swift adaption to selective pressures, whereas species with a low standing genetic variation are potentially more susceptible to rapid environmental changes (Thompson et al., 2019). The possible replacement of the large, nutritionally rich T. libellula with the smaller T. abyssorum and T. compressa is likely to negatively impact predators at the higher trophic levels. Arctic species that specialize in feeding on T. libellula, for example the little auk and polar cod, will have to adapt to a new prey spectrum consisting of smaller and less energy-rich boreal species such as T. abyssorum (Dalpadado et al., 2012;Kraft et al., 2013).
In conclusion, this study shows for the first time, contrasting molecular diversity between two congeneric species of pelagic amphipods. It contributes to a better understanding of the evolutionary processes driving molecular diversity and the adaptive potential of two ecologically important species in a changing Arctic. These results emphasize the need for further analysis of the molecular biogeography and phylogeny of key zooplankton species in the Arctic Ocean. A better understanding of population connectivity, physiology, adaptive potential and trophic ecology is crucial for formulating accurate predictions of how future zooplankton communities, species interactions and food-web structure will materialize in the Arctic and its marginal seas as a result of climate change.

ACK N OWLED G M ENTS
This study has been conducted in the framework of the Helmholtz Access funding enabled and organized by Projekt DEAL.

CO N FLI C T O F I NTER E S T S TATEM ENT
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data and sampling locations of all analysed specimens are publicly available on the BOLD depository under the project name 'ARCTH' and made available on GenBank (Accession numbers: OR210933-OR211370).